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^ ■ We review existing methods for generating long streams of 1//" noise (0 < a < 2) focusing on the digital 

filtering of white noise. We detail the formalism to conceive an efficient 1//" random number generator 
(white outside some bounds) in order to generate very long streams of noise without an exhaustive computer 
memory load. For a = 2 it is shown why the process is equivalent to a random-walk and can be obtained 
simply by a first order filtering of white noise. As soon as a < 2 the problem becomes non linear and we 
show why the exact digital filtering method becomes inefficient. Instead, we work out the formalism of 
using several 1//^ filters spaced logarithmically, to approximate the spectrum at the percent level. Finally, 
_ from work on logistic maps, we give hints on how to design generators with a > 2. The corresponding 

' software is available from http : //planck. lal . in2p3 . f r/giki/pmwiki .ph p/Sof ts/AbsRaiid] 
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Many physical systems exhibits some 1//" noise, i.e. a stochastic process with a spectral density 
having a power exponent < a < 2. Among the many fascinating questions related to it an 
apparently simple problem is how to build random number generators (RNG) with such long-range 
correlations. The classical method (see e.g. Ij) is to generate a vector of white noise (as produced 
by a standard RNG) with an appropriate correlation matrix. However, one needs to take the 
"square root" of this matrix (as with a Cholesky decomposition [2]), a method rapidly limited by 
the size of the requested sample. One must then turn on to more specialized techniques, which can 
be categorized into: 

• Digital signal processing: one generate a vector of white noise and fast-convolve it in Fourier 
space with the required 1//" shape (see e.g. One may use also a wavelet decomposition 
( [3]) which is quite natural due to the scaling properties of the process (see next part). With 
these methods samples of the order of o(lO^) can be efficiently obtained. For samples with 
longer correlations, one is limited by the memory of the computer. 



Fractal techniques: in non-linear science, 1//" noise is referred to as "fractional brownian 
noise" . The peculiar shape of the spectrum (with a power exponent) implies that the process 
is " scale- invariant" or more precisely "self-similar" (see e.g. ^^), meaning that the stochastic 
process X{t) is statistically similar (in the sense of having the same statistical moments in 
the infinite limit) to ^^^S^ for any dilatation r, where H is the Hurst exponent related to the 
1//" slope by a = 2H + 1. Some very efficient method to generate this process [B] is to shoot 



^a large collection of references is available online: 
|http : //www. nslij -genetics . org/wli/lf noise/ index . html | 



Gaussian numbers for both ends of an interval, take the average at the mean and add a new 
shot with a re-scaled amplitude, and repeat the process on each subintervals. This however 
requires knowing a-priori the limits of your sample and is also limited by the computer 
memory for very large samples. 

In some cases, these methods are insufficient, as in some modern astrophysics experiments. 
For instance, this work has been perform in the framework of the simulation of the Planck ESA 
mission [7] which is a satellite experiment, planned to be launched in 2008, whose goal is to measure 
the Cosmic Microwave Background anisotropics with unprecedent accuracy. It is composed of two 
instruments: High Frequency (HFI) and Low Frequency (LFI), the former using low temperature 
bolometers (O.IK) and the latter radiometers. At the level of precision required (~ 10~^K) the 
instruments are sensitive to 1//^ noise mostly from the cryogenic parts (HFI) and 1//" with 
a ~ 1.7 from the low-noise electronics (LFI). Given the sampling rate of the instruments (200 Hz) 
and its long duration (two years), very long streams of noise need to be generated to prepare the 
analyzes . The techniques described previously fail. 

In the following we will show how to produce rapidly such long streams of 1//" noise for an 
arbitrarily low frequency cutoff, without requiring massive memory load. 

In the 1//^ case (section [1]) we will adapt some classical numerical filtering technique (Auto- 
Regressive Moving Average, see Appendix) that allows to build an optimal generator {i.e. only 
limited by the underlying standard white one). For a < 2 (section [5]) the problem becomes much 
more intricate, because we are moving off the linear theory and the previous method cannot be 
generalized. We will then adapt a method from electronics, to approximate the problem to a very 
good precision by a set of 1/ generators and use the previous case. 

We will then give hints on how to build generators with a > 2 (section [3]) using methods based 
on logistic maps. 

The Appendix recalls some properties of the z transform which will be our main tool. 

The algorithms described hereafter are now full part of the Planck simulation programs p| and 
streams of 1 year data of noise (~ 6.3 10^) are generated in about 20 minutes on a single standard 
2 GHz processor. 

1. 1//^ noise 

1.1. Why pure 1//^ noise is random walk. 

Let us start with the pure 1//^ noise case. 

The idea is to start with a standard Gaussian generator and numerically filter the shots (called 
Xi, realizations of an x{t) process) to obtain the proper power spectrum which is simply [^: 

Syico) = \Hiju;)\^SAu;) (1) 

Since x is white S'j; = 1 , so all we need is to design a filter, here of the form: 

= ^ (2) 

therefore: 

Hij.) = -L (3) 

This is simply a first order integrator and has an equivalent discrete z transform of the form (see 
Appendix): 

Hi^-') = (4) 

Therefore the filtered signal is: 

Y{z-')^H{z-')X{z-') = f^ (5) 

^Through the article we will use the electronic notation for = —1 and work in the frequency domain uj. In 
addition we will not mention the normalization factor a of the input Gaussian generator that can be recovered 
straightforwardly by multiplying the output by 
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which according to the summation property of the z transform (Appendix) is obtained simply from: 

fe 

Vk^^Xi (6) 

1=0 

This equation represents an MA (moving average) filter; the system is built by "remembering" at 
each step its previous value (random walk) and has therefore only one state variable. It has an 
"infinite memory". 

It is slightly more convenient to write this filter in the following way: 

(i-z-i)y(z-i) = x(z-i) (7) 

which accordingly to the offset properties of the z transform is obtained in the time domain through: 

Vk - Vk-i = Xk (8) 
which represents a one state AR (Auto Regressive) filter. 

1.2. Whitening below some minimum frequency fmin 

While infinitely long streams may be theoretically interesting (does pure 1//" noise has an intrinsic 
low frequency cutoff? Is it a stationary process?...) we are concerned here with designing a 
generator for real life experiments which have finite durations. Therefore it is assumed that there 
exists some frequency /min below which the noise becomes white (which can be as long as the 
time-scale of the experiment). 

To get a white spectrum below some frequency ivq — 27r/min, we are seeking for a transfer 
function of the type: 

mM\' ^ (9) 

H{ju;) = (10) 

which represents a first order low-pass filter. 
Its inverse Fourier transform is: 

h{t) = woe"""*w(t) (11) 
{u{t) being the unit step), which sampled at a rate 1/T has the z transform: 

Hiz-^) = Y: hmz-^ - l-e^lr,-, (12) 

i 

Therefore: 

y(z-i)(l-e-"°^z-i) =cjo^(z-') (13) 

whose inverse transform is: 

Uk =uJo Xk+ e^^^^T/fc-i (14) 

For each shot Xk one just needs to keep the previous filtered value yk-i to compute the new 
one yk- 

1.3. Whitening above some knee frequency fknee 

In most experimental cases, noise becomes white also above some "knee" frequency /knee and we 
incorporate it in the following. 

For a filter with a white spectral behavior both below /min and above /knee, we symmetrize 
Eq.©: 

mju.)f . ^ (15) 

H{juj) ^ —— (16) 

JUJ + UJo 



In order to treat the numerator and denominator symmetrically, we use the bilinear transform (see 
Appendix): 

relating it to frequency via: 

U W 2w/samplG (18) 

to get the z transfer function in the form: 

//(.-) ^ ^4±^ (19) 



with: 



1 - hiz- 



l + ri 
1 + 



ai = -^—^ (20) 



1 + ro 
1 + ro 



where the reduced frequencies are ri = 7r/kncc/ /sample and tq = vr/,„i„//sampio- 
It is now straightforward to build the numerical filter from: 

r(z-i)(l - 6iz-i) = (ao + aiz-^)X{z-^) (21) 

which for any k reads: 

Uk = aoXk + aiXk-i + hyk-i (22) 

This is the equation for a simple ARMA filter: for any white shot Xk one just needs to keep 
the previous one Xk-i and the previous filtered value yk-i, to compute the new one. 

This simple filter has been implemented in C++ and run for 10^ samples. Results are depicted 
on the upper plots of Figure [T] and show that it behaves correctly. 

This method is optimal, in the sense that CPU-time is just limited by the speed of your white 
number generator. It is insensitive to whatever your maximal correlation length (/min) is. And 
obviously it does not require any memory load. 

2. 1//" noise 

2.1. Failure of ARMA models 

ARMA models were so successful in the 1//^ noise case that we are tempted to generalize them to 
the 1//" case. This is however far from trivial, mainly because we are moving off linear theories 
for which the z transform has been tailored. Let us see what happens in more details. 
For a pure 1//" noise case,we are seeking for a transfer function of the form: 

\Hiju;)\' = ^ (23) 
Its generalized ^ inverse Fourier transform is: 

and its associated one-sided z transform is: 

H{z-')^Y.h{^T)z-^ = (26) 



CD 

> 
CD 




5000 10000 
sample number 




'■5 10"^ 10"^ 10'^ 10"^ 
frequency (Hz) 



10" 



CD 

CC 
> 

CD 




5000 10000 
sample number 



N 
X 



1000.0 



100.0 



V) 

d 

CD 
T3 

o 

CD 
Q. 

CO 



10.0 




frequency (Hz) 



Fig 1. Results obtained with the generators. Upper plots: random walk (a = 2). Lower plots: 1// noise (o = 1). 
10^ samples are shot in each case. A slice of the signal in real space is shown on the left, and the power spectral 
density on the right (estimated using a classical Welsh periodogram method [2]). The following bounds were used 
for whitening: /min = 10~* Hz, f knee = 10~^ Hz. The logarithmic scale allows to check by eye the validity of the 
slope. 

Unfortunately the sum that appears on the r.h.s cannot be factorized, so that the MA filter 
Y{z^^) = H{z^^)X{z^^) leads to an explicit convolution: 



/2; 



Xk- 



(27) 



i=0 



In direct space, one can therefore produce I// noise by generating type transients at random 
initial times [TFlfTg] , 

There are however two problems with this expression: 

1. the first term z = is infinite 

2. we must keep all the past terms (xfc_i)which is rapidly inefficient even with a fast convolution 
method. 

The first problem is the reflect of the high frequency divergence of I//" noise [THj. Kasdin [TU] 
proposes an elegant way to circumvent it, by using a z transfer function of the form : 

1 



Hiz-^) = 



(1 -Z-l)"/2 



(28) 



(compare to ([S])). It allows us to build an AR filter: Y{z ^)(1 — z ^) "/^ = x{z ^) by the power 
expansion (1 — z"^)""/^ = ikz~^: 

k 

'^ikyk-i = Xk (29) 

i=0 

The ^fc tend rapidly to fc^('^+"/2) and the first term = 1 is now well defined. 

This however docs not solve the long range dependency which just reflects the fact that a 1//" 
(a < 2) spectrum has an infinite number oj state variables and that there exits no (linear) difference 
equations to represent it. 

One may think that this bad behavior is to due to a missing low frequency cutoff (o^o)- We can 
modify the transfer function (|23|) to see what happens: 

\H{J^)? = ( (30) 

^(J^) = ?^'w/2 (31) 
then using Eq. 3.382-7 [llj , the inverse Fourier transform is: 



By comparing to the no-cutoff case Eq. lfzo]) . we see that the effect of introducing cjq is to attenuate 
the long range correlation by the exponential factor: t,k ~> CfcC Still, one needs to keep a few 

times l/(ci;or) past samples which is too inefficient for a low frequency cutoff. 

2.2. Recursive 1//^ filtering 

We then turn on to an approximate method to produce efficiently long range correlated samples. 
An old well-known way in electronics (12j is to approximate the 1// spectrum by a sum of 1//^, i.e. 
relaxation processes, equally spaced on a logarithmic frequency grid, which is called the "infinite 
RC model" [TB]. Quite surprisingly little filters are needed to produce an 1//" (0 < a < 2) 
spectrum at the percent level. This method is used for instance in pink audio noise generation |14j , 
in commercial libraries [E], and even implemented on a DSP [T^ . 

We work out here a slight variant of the RC model due to Keshner [17] which allows to incor- 
porate whitening both below /min and above /knoc, by filtering one single white shot. 

We are seeking for a transfer function of the type: 

The approximate transfer function is of the type: 

Nf-l 2 2 

i^(^-)i^-n^ (34) 

i=o ^ 

where Nf is the number of filters used, and pi, Zi are the poles and zeros locations. 

Since we already designed these 1//^ filters in section [TT3l the only remaining point is to get 
their poles and zeros positions. 

From Figure [2] it is clear that the poles must be located regularly on a logarithmic gri 

logWfe-logWo 

Ap = (35) 

logpi+i = logpi + Ap (36) 



^we use a base 10 logarithm. 
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Fig 2. A band-limited 1// spectrum (full line) can be approximated by a set of 1//^ spectra properly located 
(dashed lines) . The crosses and diamonds represent respectively the position of the poles and zeros in Laplace 
space. Here we use only one filter per decade; the approximation is already good. 

in order to have a symmetric errors, the zeros must be placed at [20j : 

log Z, = \0gp^ + (3'^) 

and for continuity, from geometrical considerations, the first pole must be put at: 

f a 

logpo = logcJo + ^(l - (38) 

How many 1//^ filters do we need? While in principle, as few as one per decade is enough to 
approximate the slope to better than 5% [17], since we also wish here to reproduce neatly the /min 
and /knee shoulders, we used 1.5 filters per decade, which leads to an approximation at the percent 
level for any < a < 2 value throughout the spectrum 0. 

We implemented this algorithm, using the previous code for 1//^ noise: a Gaussian white 
number is shot and passed recursively through the set of 1//^ filters: the coefficients of Eq 
are automatically recomputed by the class for a given pole/zero and the output of each filter is the 
input to the next one. 

We illustrate on the lower plots of Figure [T] results obtained by generating 10^ samples of an 
1// noise. It required using only six 1/ filters and the increase in CPU-time with respect to a 
pure white Gaussian shot, is only of 25%: the limitation is mainly due to the white time shot, not 
the filtering stages. 

The strong feature of this method is that the number of filters used increases only logarithmically 
down to /min- The drawback is that the generator has to be "warned up" to avoid the transient 
response, on a timescale of the order of the slowest RC filter (approximately l/Zmin)- 

3. a > 2? 

While rarely encountered in experimental data, the curious reader may wonder how to generate 
an 1//" noise with a slope a > 2. From geometric arguments, the methods presented before fails. 
An interesting approach is to use the work on logistic maps. 
Consider the straightforward (non-linear) recurrence [21) : 

Xk+i — Xk + xj. (mod 1) (39) 



*The worst case being the 1// spectrum. 



where < < 1. 

This is a pure 1// noise generator (actually l//(log/)^)! 

By changing the power exponent, group theory computations show ^22J that any factor a > 2 
can be reached. By adding a "shift from tangency" e in the recurrence: 

Xk+i = (1 + e)xk + (1 - e)4 (mod 1) (40) 

one introduces a low frequency cutoff fmin- Whitening above /knee can be obtained by shooting a 
properly weighted extra white shot. 

This chaotic process is however far from Gaussian (it is an "intermittent" one, representing 
long phases of low values, followed by sudden "bursts" of high ones ) but, from the Central Limit 
Theorem, Gaussianity can be recovered by summing up several such independent generators, a 
method particularly efficient on parallel computers, but beyond the scope of this article. 

For completeness, note an interesting connection between this approach and a more empirical 
one consisting in re-scaling some Brownian motion (but in which case only < a < 2 range can 
be reached) 

Conclusion 

We have worked out the formalism of filtering white noise efficiently to produce " infinite" streams 
of band-limited 1//" noise. 

• For noise the method is optimal, i.e. is equivalent to generating standard Gaussian 
white noise. 

• For < a < 2 we use an excellent approximation which only requires about one and halfl//^ 
filters/decade, allowing the very fast generation of 1//" noise over an arbitrarily large fre- 
quency range. 

The corresponding (C+-I-) software can be downloaded from: 
[htt p : //planck. lal . in2p3 . f r/wiki/pmwiki .php/Sof ts/AbsRand| 

We have focussed mainly on this method because it is lightweight and very efficient for very long 
time streams. For more general applications, before choosing a dedicated algorithm ask yourself 
the following questions: 

1. Do I want to generate a Gaussian noise? 

2. What a range do I need? 

3. Do I have a low frequency cutoff /min below which noise is white? Note that the upper /knee 
frequency can always be added by shooting an extra white noise (properly weighted). 

4. On how many decades do I need the logarithmic slope? 

5. What is my hardware support (mainly memory load)? 

Table [T] then gives you some hints about which algorithm you can use. 



Method 


Gaussian? 


a range 


fmin 


CPU/load 


ref. 


FFT/FWT 


yes 


any 


yes 


moderate 


e.g. [3]/ [3 


Random pulses 


yes 


< a < 2 


yes 


high 


[H], [IS] 


ARMA filtering 


yes 


< a < 2 


yes 


high 


[TD], section O 


Sum of filters 


yes 


< a < 2 


yes 


low 


[H], [g, [H], section O 


Random midpoint displacement 


yes 


1 < a < 3 


no 


low 




Brownian scaled motion 


no 


< a < 2 


no 


low 


m 


Logistic maps 


no 


any 


yes 


low 





Table 1. Summary of today's existing mctliods to generate discrete l/f"' noise. The column "/min " refers to the 
(possible) existence of a frequency limit in the method, beyond which noise is white. "CPU/load" is indicative (it 
should be considered as relative to the different methods) and only significant for very long correlation lengths. 



Appendix 

In this appendix we recall a few results of the theory of automatism that are relevant to this 
analysis. The reader is referred to any book on linear theories (as [23) for more details. 

The z transform is a tool similar to the Fourier transform, but for the discrete domain. To a 
causal series of values {xi=o,i,..} we associate its z transform by: 

X{z-'^) = Xo+ XiZ-^ + X2Z^'^ + ... (.1) 

where z is a complex number. 

The following properties can then be demonstrated: 

linearity Ai{xJ + AajyJ ^ XiX{z-^) + A2r(z-i) 

offset {x^-k} ^ z-''X{z-^) (.2) 

summation {EU^.}-^^ 
convolution {YX=o ^kVi-k} ^ X{z-^)Y{z-'^) 
A time signal x{t) sampled at a period T is described as: 

x*{t) = ^x{iT)5{t-iT) (.3) 
and we can associate to the samples Xi = x{iT) a z transform. Some usual transforms are: 

Sit) ^ 1 (.4) 
u{t) ^ ^.- = __^ (.5) 

OO ^ 



1=0 



An interesting property is related to energy conservation. The power of the sampled signal is 
related to its z transform: 

|X*(jw)|2=X(z-i)X(z)|,^,,.. (.7) 

The link between \H*{jLu)\ and |_ff(ja;)| being not convenient, it is sometimes interesting to use 
the w bilinear (or Tustin) transform which is defined as 

1 — 2 , , 

1 + z 

It can be shown that it is related to the frequency of the signal oj by 

w = tan(wr/2) (.9) 



In the limit of fast sampling loT — > (which is most often the case since we work below Nyquist 
frequency) the w transform can be identified with oj, giving a convenient way to determine the z 
transform coefficients directly from the signal transfer function. 

Finally, from the properties (f72|) , a rational polynomial z transfer function of the form 

leads to the following filtering in the discrete time domain: 

Vk + biUk-i + ... = a^Xk + aiXk-i + ... (.11) 



when only oq is different from 0, it is an called Auto regressive (AR) filter (also know as Infinite 
Impulse Response). If all the &i's are null, this is a Moving Average (MA) one (also know as Finite 
Impulse Response). In the general case this is called an ARMA filter. 
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